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INTRODUCTION 



I . 

Recent research at the NPS in atmospheric global predic- 
tion has been primarily concerned with the five-level global 
primitive equation model initially designed and programmed 
by Dr. F.J. Winninghoff and further developed by Elias (1973), 
Mihok (1974), McCollough (1974) and Maher (1974). Their 
numerical experiments proved valuable but pointed out the 
inflexibility of the model for research purposes. 

Arakawa and Mintz (1974) described an atmospheric global 
prediction model with a new finite difference formulation 
and greater vertical resolution. The model employs a new 
horizontal distribution of variables to improve the geo- 
strophic adjustment process. The purpose of this research 
was to program and check out a modified version of the new 
Arakawa scheme. The fundamental approach was to develop a 
flexible program for use in a number of projects including 
modeling of tropical circulation. Eventually the model's 
performance can be compared with the performance of FNWC's 
global model. 

The model was initially developed with two levels and a 
variable grid to aid in debugging. The atmosphere was 
treated as adiabatic and frictionless and the water vapor 
continuity equation was not considered. The initial condi- 
tions for all experiments were analytically derived in the 
same manner as Maher (1974). The advantages of analytic 
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initial conditions were significant in reducing computer time 
during check out of the model. 
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II. MODEL DESCRIPTION 



The differential equations are essentially the same as 
those described by Arakawa (1972). The integration scheme 
was carried out on a staggered, spherical, sigma coordinate 
system. For an adiabatic, frictionless model the primitive 
equations form a closed set. No sink or source terms were 
considered in the course of development or check out. 



A. VERTICAL COORDINATE SYSTEM 

The model uses the non-dimensional sigma coordinate sys- 
tem as described in Haltiner (1971). The two levels divide 
the troposphere in half and the tropopause is assumed iso- 
baric. The sigma coordinate is defined as 




where p is the pressure, p^ the constant tropopause pressure 
and 7 t is the terrain pressure. The terrain pressure is fur- 
ther defined as 

v E p s " P t ' (2.2) 

where p is the surface pressure. It follows from equation 
s 

(2.1) that 

a = 0 at p = p t , 

(2.3) 

a = 1 at p = p s , 

which are the vertical boundaries of the coordinate system. 
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Figure 1. The sigma (o) coordinate system as used in the 
model. is 200 mb. 
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The boundary condition at a=l is a=0 and at a= 0 it is as- 

• • 

sumed that cr=0 where o is defined as dcr/dt. 



B. PRIMITIVE EQUATIONS 

The continuous form of the primitive equations is de- 
scribed in the orthogonal curvilinear coordinates E, and n 
where E, = X (longitude) and n = <p (latitude). An area ele- 
ment is — A£An and the lengths of the sides are — and — , 
mn m n 

where 



1 , , 1 
— = a cos d> and — = a. 
m n 

The difference element is illustrated below. 




m 



The mass continuity equation is 



3 , it , 9 , u. ,9 . v. ,9 jo. 
9t^mn^ 9£^ U n^ 9a^mn^ 



= 0. 



(2.4) 
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The ^-component of the horizontal equation of motion is 



8 , tt . , 8 /iru, v . 8 , ttv v , 8 ,t\ a N 

8t^mn u ^ 8C n u ^ 8n m u ^ 8a^mn u ^ 



, f , , 8 1 8 1., , 7T r 8$ , 8 7T IT -r, /o CN 

- t — + (v-r-? — - u-z )]ttv + — [-r— + aa -r? = — F_ . (2.5) 

mn 3^ n 8n m n 8£ 8C mn £ 

The n-component of the horizontal equation of motion is 



— (_!L v) + + i-(IXv) + 

8t^mn ' 8^ n ; *r\ K m ' 8a^mn v; 



r f , / 3 1 8 L, , tt r 8 <J> 8 tt. 

+ [ + (v-rrr — - Ut: ) ] 7 TU + — [-r— + a a tt— ] 

mn 8£ n 8qm m 8n 8n 



— F 
mn n 



( 2 . 6 ) 



The first law of thermodynamics is written in the following 
manner : 

c p| = «« ( «. < 2 - 7 > 

where 

a) = -^2 = ira + a(-|r- + V*V)ir (2.8) 

Qt d t 

and the flux form is 

I^TTCpT) + V a -(TTV C p T) + |_ (7r aC p T) = ir(wa + Q). (2.9) 

Using the relationship for potential temperature, 

T = e(^-) K , (2.10) 

Po 

the final form in curvilinear coordinates is 
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( 2 . 11 ) 



The equation of state is 



a 



RT 

P 



( 2 . 12 ) 



and the hydrostatic relationship is 



6$ = -ira6a . 



(2.13) 



A complete set of symbols for the above equation may be 
found in the front of this report. 

C. VERTICAL INDEX 

The index k is used to identify the vertical levels. At 

the upper boundary k=0, p=p and at the lower boundary 

k=K+l and p=p . The variables V and T are carried at the 
s 

odd levels while tic is carried at the even levels. The fol- 
lowing definitions are required to complete the vertical 
system: 




(2.14) 




(2.15) 



where the summation is over odd k (see figure 2). 
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VERTICAL INDEX 



Figure 2. In the above figure k is a variable vertical 
index, sigma (o) is the dimensionless vertical coordinate, 

V is the horizontal vector velocity, ir is the terrain pres- 
sure, c is the vertical velocity, T is the temperature and 
4> is the geopotential. 
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D. HORIZONTAL GRID AND DISTRIBUTION OF VARIABLES 

The selection of the horizontal distribution of variables 
is a function of two distinct processes. The first is proper 
simulation of the geostrophic adjustment process and the 
second is proper simulation of slowly changing quasi-geo- 
strophic motion after it has been established by geostrophic 
adjustment. As shown by Winninghoff (1968) geostrophic ad- 
justment depends on how the variables are distributed over 
the grid points. Winninghoff used the following equations, 
which are the simplest ones in which geostrophic adjustment 
can take place, to demonstrate five possibilities for the 
placement of the dependent variables: 



du 

dt 



fv + g 



ah 

ax 



= o , 



(2.16) 



dv 

dt 



TTr + fu + 



ah „ 

g— = 0 

oy 



(2.17) 



dh , . ,3u , 3v, „ 

dt + h <33 + sj) = 0 • 



(2.18) 



The above equations represent an incompressible, homogenous, 
non-viscous, hydrostatic, rotational fluid with a flat bot- 
tom and a free surface, where u and v are the velocity com- 
ponents, h is the depth of the fluid, f is a constant 
Coriolis parameter, t is the time, x and y are the horizon- 
tal coordinates and g is gravity. The five possible distri- 
butions of the dependent variables h, v and u are shown in 
figure 3. Scheme B was used in the Mintz-Arakawa two-level 
global circulation model (Langlois and Kwok, 1969) and 
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Figure 3. The above figure represents the placement of the 
dependent variables h, v and u where h is the depth of the 
fluid, u and v are the velocity components and d is the grid 
distance (for the simplest case of geostrophic adjustment). 
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scheme E is used in the global model now under further de- 
velopment at FNWC (Maher, 1974). 

In the one-dimensional case it was shown that both 
scheme B and C adequately simulated geostrophic adjustment. 
In the two-dimensional case, however, scheme C was shown 
most satisfactory to simulate geostrophic adjustment, except 
where X R /d is less than or close to one. The quantity X R is 
the Rossby radius of deformation and d is the grid distance 
(Arakawa and Mintz, 1974). The Rossby radius of deformation 
is 

x R = /HT/f 

where H is the mean value of h, g is the acceleration of 
gravity and f is the constant Coriolis parameter. The con- 
dition X R /d £ 1 is an abnormal case, therefore, the horizon- 
tal distribution of variables is based on scheme C. (See 
figure 4 . ) 

E. TIME DIFFERENCING • 

The time differencing is carried out in thirty minute 
sequences. The initial step in each sequence is a two part 
Matsuno scheme represented by the following notation: 

F* = F* + At •|| t (Forward), (2.19) 

O t 

* 

F t+1 = F^ + At -T 7 * (Backward) . (2.20) 
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Figure 4. The above figure represents the horizontal distri- 
bution of dependent variables with polar modification. P in 
the ij index system represents the north pole, 1 represents 
the south pole and i represents a meridian. n represents 
the variables T, $ and it carried at "iT-points." u and v are 
the horizontal components of velocity. 
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F is a vector representing the dependent variables. The 
superscript t represents the time step, the superscript * 
represents the results of an intermediate step and At is the 
time interval of a single step. The remaining time steps in 
a sequence use the leapfrog scheme as follows: 

F t+1 = F t_1 + 2At . (2.21) 

At the end of each sequence provision is made for the calcu- 
lation of the source terms. The present model doesn't in- 
corporate source term calculations. 

F. AVERGING THE PRESSURE GRADIENT AND ZONAL MASS FLUX NEAR 

THE POLES 

A problem that may arise at higher latitudes is computa- 
tional instability. Computational instability is a result 
of the convergence of the meridians to the poles causing 
greatly reduced grid distances along a latitude circle. 
Therefore, some technique must be employed to eliminate the 
stability problem. One technique involves reducing At in 
higher latitudes but it would impose serious programming 
difficulties and require more computer time. Another tech- 
nique would be to reduce the number of grid points in a 
latitude circle at higher latitudes, however, it also would 
impose programming difficulties. A third technique is zonal 
smoothing which was developed by Arakawa and Mintz (1974). 
Smoothing preserves the integrity of the horizontal grid 
while maintaining a constant time step. Zonal smoothing is 
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used in this model and also by FNWC to eliminate computa- 
tional instability at high latitudes. 

Zonal smoothing involves expanding the pressure gradient 
and zonal mass flux into a Fourier series and reducing the 
amplitude of each wave component by a factor 



S 



D cos 4> 

C o A t sin ( m^d ) 



( 2 . 22 ) 



where 

S = Stability Coefficient, 

D = Grid distance at the equator, 

Co = Phase speed of the fastest gravity wave, 

At = Time step, 

m^ = Wave number , 

d = Grid distance in degrees. 

If the value of S is greater than one, smoothing is not re- 

quired (Arakawa and Mintz, 1974). As pointed out by Arakawa 
the smoothing operation does. not smooth the fields of vari- 
ables, because it is simply a generation of multiple point 
difference quotients. The bar operators in Chapter III in- 
dicate zonal smoothing. 



G . PROGRAM FORMAT 

The program format was patterned after the Mintz-Arakawa 
two-level global model as described by Gates, et al. (1971) 
and Langlois and Kwok (1969) and modified to include scheme 
C distribution of variables and a Matsuno-leapf rog time in- 
tegration scheme. Figure 5 gives a simplified flow diagram 
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Figure 5. A simplified flow diagram for 5-day forecast 
using analytic initial conditions. 
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for a five-day forecast. The main program controls the 
overall time period of the forecast as well as input and 
output. Time integration for each 30-minute sequence is 
controlled by subroutine Timestep. The forcing functions 
are calculated in Comp 1 and Comp 2 with zonal smoothing ap- 
plied as necessary. The remaining subroutines are peripheral 
to the main flow of the model. 

A main consideration in the program was flexibility, 
therefore, the vertical structure along with the horizontal 
grid are variable. The model can be used to simulate sever- 
al different horizontal and vertical structures including 
the FNWC's 5- level global model. 
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Ill . FINITE DIFFERENCE EQUATION 



The finite difference equations given in this chapter 
were developed by Arakawa and Mintz (1974). The advection 
terms of the equation of motion were modified to eliminate 
the diagonal flux calculations. The short-term performance 
should not be affected by elimination of the diagonal flux 
terms and considerable computer time was saved. 

The notation used in this chapter represents the specif- 
ic variable centered in the ij index system. Figure 6 con- 
tains an example of "ir-centered" and "u-centered" notation. 
IT-centered notation is used for the continuity equation and 
the thermodynamic equation in which the dependent variable 
is carried at a iT-point. U-centered notation is used for 
the ^-component of the equation of motion and v-centered 
notation is used for the ri-component of the equation of 
motion. Refer to the front of the report for a complete 
list of symbols. 

A. CONTINUITY EQUATION 

The following form is used for the continuity equation 
given in equation (2.4): 



^^i i k k k k 

+ F i +i ,j - Vi.j + G i,j +i - 



+ — — r- (S k+1 - S k *) = 0 , 
Ao k 1>J 



(3.1) 
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^CENTERED 




Figure 6. An example of the notation used to describe the finite dif- 
ference equations. The continuity equation is described on a "TT-centered" 
grid in which the F and G symbols are flux calculations in their re- 
spective directions. The M u-centered" grid is an example used to de- 
scribe the ^-component of the equation of motion where F U and g U are 
also flux calculations. 
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where 



F i + 2,3 n ^i+i, j ^i+l,j + *i,j^ ’ 



(3.2) 



(bar indicates zonal smoothing) 






. (7T . . . - + IT . . ) , 

m 'i'j+x i, J+l i,J 



(3.3) 



n . = tt (ASAii) 

i,j i,j ^ ran . 

J J 



s . . = n . . a . . 
1,J 1,J 1,J 



The tendency equation is 



an. 

TT 



i , j _ 



K, 



- r ' c 

k=l 



i+i, j 



i — 2,3 



+ - G^ NArr* ( 'X A 

• U -j -i+i U -j + 

1 > «J^2 1>J t 2 



The vertical motion equation is 



s k+1 

i, J 



-J; (f ih,j - F i-i,j + 



~k u k 
G. . i)Ao 

1 > J“2 



- a 



k+l 9II i.j 

at 



(3.5) 
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B. PRESSURE GRADIENT Tkum AND HYDROSTATIC EQUATION 
1. Hydrostatic Equation 

The geopotential for a level K is given by 



,K 
$ . . 
1 , 3 



i, J 



K, 

k=l i , J 



k 



Acr 



i, J 



K e 2, k+1 ^ ~k+l r/ p k+2 K “ k K 

k=l 



C 0* i [(^ ) - (£-) ] , (3.6) 



P i, J 



i, J 



i, J 



g 

where $ is the geopotential at the surface, 



;k+l _ £n 0 k - £n 0 k+2 



0 k+2 0 k 






Given $ , the geopotential for the remaining levels is found 
using the following relationship: 



k+2 k 



k k 



<!> k . - <I> k+2 = C [(2 ) - (£-) ] e k+:l ; 

1 » J P P 0 ±> j P 0 J j 



(3.7) 



2 . q-component of the Pressure Gradient Force 

The pressure gradient force in the meridional equa- 
tion of motion as given in equation (2.6) is 

- S 4 4[(,, i i++ + "i i+i - *i i-D 

111 1,J^2 l,J- 2 l,J- 2 



+ ((™a) k . i + (Tiaa). . x )(tt . x - tt . x )] , (3.8) 

l,J + 2 1>J“2 l,J + 2 1,J _ 2 
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where 



. , ,k k RT ?, j+ 4 

( ” a) i,j + 4 "i.j+i 0 



3. ^-component of the Pressure Gradient Force 

The pressure gradient force in the zonal equation of 
motion as given in equation (2.5) is 



2 [(it. . + tt . . )($ k - . - <I> k .) 

n i+l, J i,J i+l, J i, J 

J 






(3.9) 



O MO v *l 7, ^ T, T‘T™ f THT TTVTP Q ■ 

The zonal momentum flux as given in equation (2.5) is 



-^r(II U • u k .) + i[ F U . i • ( u . , - . + u. .) k 

3t VJl i,j 1,3 L i+i,J i+l, J i,J 



^i-2,j^ U i,j + U i-1, j ^ + g i, j+l^ U i, j+1 + U i , j ^ 

- g U . l( u. . + u. . 1 ) k ] + JL. i[§ U ’ k+1 (u k+2 + u k .) 
S i, J~i V i,J i,J-l J . k i,j i,J i,J 



-§ U,k-1 (u k . + u k ~ 2 )] , 
i,J i,J i,J 1 



(3.10) 



where 






= 4(n i+i + n. i .) , 

1 + 2 ,J 1 -^>J 
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5^ . = 4(S_, . + s. , .) 

i ,3 ■ i+i,j i-£,j 



F U and g U are defined as follows: 



F U,j 5 * (F M,j + l + 2F iH,3 + F i + *.J-1> ’ 



(3.11) 



e u 

i , J + 2 



= x 



i (G . , + G A , +1 + G. i . + G x . +1 ) , (3.12) 

1 + 2,J 1-5, J 1-2, J + 1 



where 



* 

F. . 

1 , 3 



= i 
- 2 



(F i+1 4 + F ) , 

2 ) J 2 j J 



(3.13) 



I. . - 2(G. + G. . i) 

1,J 1,3+5 1 ,3-h 



(3.14) 



The meridional momentum flux is similar to equation (3.10) 
with u replaced by v. F v , g v , n v and S v take the following 
form: 



^i+4,j 4 ^ F i+l,j+i + F i,j+i + F i , j-2 + F i+1 , j-i ^ 



g i,j+i 4 ^ G i+l,j+i + 2G i,J+i + G i-l,j+l^ 



n i i E i(I1 i 1+i + n i i-*> ’ 

1 > J 1 > 2 1 > J 2 



• V „ • • 

s: . e i(s. + s. . x ) 

i,J i,J+l 1,3-4 



F is defined by equation (3.13) and G by equation 

/ 



(3.15) 

(3.16) 

(3.17) 

(3.18) 
(3.14). 
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D. CORIOLIS FORCE 



A variable C. .is defined at TT-points by 
i > J 



c k . = f. .(A££H) _ a (u . + u. i .) k 

i, J ij nrn ^ 2V i+i,j i-i,J 



t< ^ ) - 1 ■ 
m j+i m j-j 



where f. .is the Coriolis parameter. 

1 > J 

The zonal-component of the Coriolis force as given in 
equation (2.5) is 





c k 

i+i, j 


^ V i + i. j+i 


.k 

+ V i+ A i_A) 

A ' 2 » J 2 




+ "i-i.j 


C k x . 
i-i, J 




. ,k, 

+ v i_x i_i) 1 • 

1 2 } J 2 


(3.19) 



The meridional-component as given in equation (2.6) is 
-j_i C. • i (u. +i • i + u • jl .si) 

1»J 2 1>J 2 1^2 > J 2 1 2>J 2 



+ 11 i i+i C i i+A ( U i+A i+i + U i_i i+i) k ] 
1| J t 2 1>J^2 1'2>J t 2 1 2) J t 2 



(3.20) 



E. THERMODYNAMIC ENERGY EQUATION 

The thermodynamic equation corresponding to equation 
(2.11) may be written as 
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h ‘"i.j T i,j> + 



T + T 

( ■ i+1 >J i.j ) k 



- F k x . ( 



T . . + T . T + T 

i-iJ — ' i.lA ! . J ( i > J + 1 i > J 

2 ; i,j+i C 2 ; 



-G . , ( 

l ,J~i 



— i->J i + — - .[ S k+1 i J ) 

2 Aa k i'j P 0 i,J 

k-1 



• k— 1 i K ^k-1 1 k i 



+ i{(u ^ ) 1+ , ' ( " aa> ti,j +( ™ a >lj )( \ + i,j - "i.p 

X ^2 > J 



+ (uAp-) ((naa) k . + (TTaa) k - .)(ir. . - it . 1 .) 

n ,-_i -j i,J 'i-l, J i,J i-l,3 

1 2 y J 






+ (v~^) ((iraa) k . + (iTaa) k . 1 )(tt. . - ir . . .. ) } 



m . . 



l, J~i 



1 > J 



i,j-l i , j i > j - 1 ' 



+ n. . Q k .] 

i,3 i,J J 



(3.21) 



F. POLAR MODIFICATION 

The poles of a spherical coordinate system are singular 
points and polar velocity components cannot be defined. 
Terrain pressure (ir) at the poles can change as a result of 
meridional flux at all points P-| and P+i where v is carried. 
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The continuity equation is modified at the poles by omitting 
all undefined flux terms. Each pole is treated as a series 
of points in order to facilitate programming (see figure 7). 
Each index i,P represents the shaded area shown in figure 8. 
The continuity equation is integrated at each point repre- 
senting the pole and then averaged to determine a polar val- 
ue. The thermodynamic equation and the vertical velocity 
are treated in a similar manner. The advective terms in the 
equation of motion are given special treatment. The remain- 
ing terms in the equation of motion are unchanged except in 
the case where they are undefined and omitted. 

1. The Advective Terms in the ^-component of the Equa- 
tion of Motion 

The polar modification of the ^-component of the 
equation of motion is shown in figure 9. It was modified 
from the Arakawa scheme by eliminating the diagonal flux 
terms. The following form is used: 





+ u 



i,P-2 





(3.22) 



where 
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a 



i 

a 



CN 
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Figure 7- The rectangular grid used to represent the globe treats the poles as a series of points and 
requires averaging to determine the polar value of tt, T and S. The u, and v equations of 
motion are modified at the poles and treated in a separate manner. 



7 



Figure 8 



7T 




Each index i,P in the continuity equation is 
represented by the shaded area. tt at the poles 
can only change as a result of G. The thermo- 
dynamic equation and vertical velocity are treated 
in the same manner. 
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Figure 9. The polar modification of the u equation of motion- 
Z 7 * 11 and g u are flux terms. The shaded portion 
represents the area associated with each variable 
u. 
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g U given by equation (3.12) 





2 . The Advective Terms in the n-component of the Equa- 
tion of Motion 

The polar modification of the n-component is repre- 
sented in figure 10. It was modified from the Arakawa 
scheme in a similar manner as equation (3.22). The follow- 
ing form is used: 




k 





(3.23) 



where 





g V given by equation (3.16) 





42 




Figure 10. The polar modification of the v equation of 

* v * V 

motion. F and g are flux terms. The shaded 
portion represents the area associated with each 
variable v. 
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IV. ANALYTIC INITIALIZATION 



The initial conditions were similar to those used by 
Maher (1974). Analytic initial conditions offer several ad- 
vantages in the initial evaluation of a model; they simplify 
the otherwise difficult task of balancing and interpolating 
initial conditions from constant pressure surfaces to sigma 
(a) surfaces, they allow the angular phase speed to be esti- 
mated from a non-divergent model and they allow the simula- 
tion of atmospheric states that otherwise would be difficult 
to simulate. 

A. ANALYTIC BALANCING 

The initial terrain pressure and velocity fields were 
obtained from the stream function solution to the linearized, 
non-divergent vorticity equation (Haurwitz, 1940). The com- 
plete vorticity equation was later solved by Neamtan (1946). 
The stream function, ip , may be written as follows: 

m k 

ip = A sin(m^X - vt) sin <J> cos <j> - Ba 2 sin 4> , (4.1) 

where A and B are constants, v is the angular velocity, m^ 
is the wave number and a is the radius of the earth. The 
angular phase speed is given by 

_v_ = N(N+1) - 2 _ 2ft (4 2) 

m k C N(N+1) N(N+1) ’ V ' 

where v/m^ is the angular phase speed, N=m k +1 and ft is the 
angular velocity of the earth. Harmonic waves defined by 
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the stream function will move with a constant angular veloc- 
ity without changing shape assuming a baratropic, non-diver- 
gerit atmosphere (Haurwitz, 1940). However, the equations 
used in this model are for a divergent atmosphere and will 
give a smaller contribution from the second term in the equa- 
tion (4.2), especially with small wave numbers. 

Equation (4.1) is used as the forcing function to obtain 
the geopotential from the non-linear balance equation 



(Phillips, 1959). The geopotential perturbation, $ , may be 
written 



$ = a 2 A(<{>) + a 2 B(<j)) sin m^A + a 2 C(c{>)(2 sin 2 m^X - 1) (4.3) 



where 




[(m k +l) cos 2 <j> + (2mJ - m R - 2) - , 



2m 2 

Xr 



(4.4) 



(m k +l)(m k +2) 



2(fi+B)^ 



m 



k <j> [(m 2 + 2m k + 2) 



(m k +l) 2 cos 2 ] , 



(4.5) 




(4.6) 
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B. ANALYTIC WINDS 



The wind components were obtained as follows: 



u = — I M 

a a<f> ’ 



(4.7) 



v = 






a cos <J> 9A. ' 



(4.8) 



Performing the operations indicated above the u and v com- 
ponents may be written 



1 m k +1 

u = - —[A sin(m A - vt ) cos <f> - m A 

a k k 



m -1 

sin(m, X - vt) cos <J> sin 2 <f> - Ba 2 cos 0 ] , 



(4.9) 



v = —[A m, 
a k 



o t ^ A 

u-u y 



m k _1 

COS d) COS (in. X 

k 



. .-u \ 1 

“ V/ u y J 



(4.10) 



C . TEMPERATURE 



The temperature field was derived in accordance with the 
NACA standard atmosphere as follows (Haltiner and Martin, 
1957) : 

T(°K) = 288 - 0.0065Z, (4.11) 

Z (meters) = 44308 [1 - ( ^ 3 ^ 5 )° 19 ° 23 ] , (4.12) 

where Z < 10769 meters. 
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V. MODEL PERFORMANCE 



The evaluation of the model was a two-step process. The 
first step involved the use of a coarse grid and was essen- 
tially a scheme to minimize computer time while debugging. 

The second step involved the use of a finer grid for better 
resolution. The majority of the experiments presented in 
this chapter were run with the fine grid. 

The coarse grid consisted of 16 points N-S which gave a 
12 degree separation between grid points. Since strictly 
analytic initial conditions were used it was only necessary 
to integrate over one wave using cyclic continuity and there- 
fore a 10-point E-W grid was adequate. The E-W grid dis- 
tance was a function of wave number. Wave numbers 4, 8 and 
12 were used in the experiments. Wave number 4 gave an E-W 
grid spacing of 9 degrees, wave number 8 gave a grid spacing 
of 4.5 degrees and wave number 12 gave a grid spacing of 3 
degrees. The fine grid used 46-points N-S which gave a grid 
spacing of 4 degrees. The E-W grid was similar to the 
coarse grid. 

In all cases the model was run with two levels, a flat 
earth and no source or sink terms. The time step (At) was 
six minutes. Fourier analysis of the surface pressure field 
was used to compute phase speed and wave amplitude in all 
but experiment III in which a Fourier analysis of the wind 
field was computed. Calcomp charts were plotted at 12 hour 



47 



intervals and consisted of one wave of the surface pressure 
field over the complete N-S grid. The phase speed and wave 
amplitude were initialized through the constants A and B 
given in equation (4.1), where A is the amplitude of the dis- 
turbance and B is the amplitude of the mean flow. 

Experiment I . Wave numbers 4, 8 and 12 were used to de- 
termine the phase speed over the coarse, 16-point grid. The 
phase speed was set at 10 degrees per day by adjusting the 
amplitude of the mean flow (B). The graphs of phase angle 
(degrees longitude) vs latitude are shown in figures 11, 12 
and 13. The observed phase speed approached 20 degrees over 
48 hours for all three cases. 

Experiment II . This and all subsequent experiments were 
run over the 46-point N-S grid. External gravity waves w p re 
simulated with wave number 4 by setting the mean flow, the 
u and v components of the wind and the Coriolis terms to 
zero. The theoretical period of the gravity wave was calcu- 
lated using the following approximation: 



T = 2(3.14159) 
P " C K p 



(5.1) 



where 



K = ^N(N+1) 
P a 



C is the zonal phase speed (set at 300 m/sec), is the 
two-dimensional wave number, N=m^+1, a is the radius of the 
earth and T^ is the period. The above approximation gives a 
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period on the order of 7 hours for wave number 4. The ob- 
served oscillations of the surface pressure field as a func- 
tion of time are given in figure 14. The period is on the 
order of 6-7 hours. 

Experiment III . In this experiment pure advection was 
simulated with wave number 4 by removing the Coriolis terms, 
the pressure gradient terms and the vertical advection 
terms. The surface pressure and the temperature were kept 
constant. The u and v components of the wind were Fourier 
analyzed to determine phase angle and amplitude. The theo- 
retical advection was determined by taking the mean zonal 
component of the wind (u) over a given time period. A per- 
centage of the actual advect ion/theoretical advection vs 
latitude for a 6-hour period is shown in figure 15. 

Experiment IV . This experiment was designed to analyze 
the 3 effect as given in equation (4.2). It can be shown 
that the 3 effect is represented by the last term in equa- 
tion (4.2). The mean flow, B in equation (4.2), was set to 
zero which eliminates the large N-S pressure variation. The 
3 effect is most pronounced for low wave numbers. Wave num- 
bers 4, 8 and 12 were run. According to equation (4.2) with 
B set to zero, wave number 4 should retrogress at 24 degrees 
per day, wave number 8 should retrogress at 8 degrees per ■ 
day and wave number 12 should retrogress at 4 degrees per 
day. The observed retrogressions are shown in figures 16, 

17 and 18. The actual retrogressions were somewhat less 
than expected due to the mean divergence which reduces the 
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magnitude of f? . The analyzed and 48-hour forecast surface 
pressure fields are shown in charts A through F. 

Experiment V . Wave numbers 4, 8 and 12 were run with 
the mean flow adjusted to give a phase speed of 10 degrees 
per day. As the wave number increases the second term in 
equation (4.2) has less contribution and the phase speed ap- 
proaches the initialized phase speed. Figures 19, 21 and 23 
show the observed phase angle vs latitude for wave numbers 
4, 8 and 12 respectively. As expected, the phase speed in- 
creased as wave number increased. The wave amplitude vs 
latitude are shown in figures 20, 22 and 24. In all three 
cases the wave form remained essentially constant over the 
48-hour forecast period. The analyzed and forecast surface 
pressure fields are shown in charts G through L. 
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LATITUDE 




Figure 11. Phase angle (degrees longitude) vs latitude, for 
16-point N-S grid, wave number 4, phase speed 
10°/day and A = 7.0 x 10 7 . (Latitudes with zero 
wave amplitude are not included and time is 
given in hours. ) 
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LATITUDE 




Figure 12. Phase angle (degrees longitude) vs latitude for 
16-point N-S grid, wave number 8, phase speed 
10°/day and A = 1.6 x 10 8 . (Latitudes with zero 
wave amplitude are not included and time is 
given in hours.) 
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L AT I T UDE 




Figure 13. Phase angle (degrees longitude) vs latitude for 
16-point N-S grid, wave number 12, phase speed 
10° /day and A = 7.0 x 10 7 . (Latitudes with zero 
wave amplitude are not included and time is 
given in hours.) 
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Figure 14. Terrain pressure amplitude vs time (hrs) for gravity wave oscillations, 
wave number 4, B (mean flow) = 0, u and v components of the wind = 0, 
Coriolis = 0 and A = 7.0 x 10 7 . 




Figure 15, Percentage of actual advect ion/theoretical advection vs 

latitude for wave number 4 with A = 7.0 x 10 7 , phase speed 
10°/day, Coriolis = 0, pressure gradient = 0, vertical 
advection = 0 and terrain pressure and temperature constant. 
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Figure 16. Phase angle (degrees longitude) vs latitude for 46-point N-S 
grid, wave number 4, phase speed -24°/day, B=0 (no mean 
flow) and A = 7.0 x 10 7 . (Latitudes with zero amplitude are 
not included and time is given in hours.) 
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Chart A. 




«.90*S 

Initial surface pressure analysis for 46-point N-S grid, wave 
number 4, phase speed -24°/day, B=0 (no mean flow) and 
A = 7. 0 x 10 7 . 
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number 4, phase speed -24° /day, B=0 (no mean flow) and 
A = 7.0 x IQ 7 . 
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PHASE ANGLE (DEGREES LONGITUDE) 

Figure 17, Phase angle (degrees longitude) vs latitude for 46-point N-S 
grid, wave number 8, phase speed -8° /day, B=0 (no mean flow) 
and A = 1.6 x 10 8 . (Latitudes with zero amplitude are not 
included and time is given in hours.) 
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~90*N 




-90*S 

Chart C. Initial surface pressure analysis for 46-point N-S grid, 

wave number 8, phase speed -8° /day, B=0 (no mean flow) and 
A = 1.6 x io 8 . 
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~90*N 




90*S 

Chart D. 48-hour surface pressure forecast for 46-point N-S grid, 

wave number 8, phase speed -8° /day, B=0 (no mean flow) and 
A = 1.6 x 10 8 . 



61 



70 




ui 

Q 

D 



- 10 ° -5 0 

PHASE ANG L E ( DEG RE ES LONGITUDE) 



Figure 18. Phase angle (degrees longitude) vs latitude for 46-point N-S 
grid, wave number 12, phase speed -4°/day, B=0 (no mean flow) 
and A = 7.0 x 10 7 . (Latitudes with zero amplitude are not 
included and time is given in hours.) 
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wave number 12, phase speed -4° /day, B=0 (no mean flow) and 
A = 7.0 x IQ 7 . 
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~90*N 




_90*S 

Chart F. 48-hour surface pressure forecast for 46-point N-S grid, 

wave number 12, phase speed -4° /day, B=0 (no mean flow) and 
A = 7.0 x IQ 7 . 
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0 ° 5 ° 10 ° 15 ° 20 ° 

PHASE ANGLE (DEGREES LONGITUDE) 

Figure 19. Phase angle (degrees longitude) vs latitude for 46-point N-S 
grid, wave number 4, phase speed 10° /day and A = 7.0 x 10 7 . 
(Latitudes with zero amplitude are not included and time is 
given in hours.) 
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Figure 20. Terrain pressure amplitude vs latitude for initial field 
and 48-hour forecast, wave number 4, phase speed 10° /day 
and A = 7.0 x 10 7 . 
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90* S 

Chart G. Initial surface pressure analysis for 46-point N-S grid, wave 
number 4, phase speed 10° /day and A = 7.0 x 10 7 . 



67 



90*M 





Chart H. 48-hour surface pressure forecast for 46-point N-S grid, 
number 4, phase speed 10°/day and A = 7.0 x 10 7 . 
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Figure 21. Phase angle (degrees longitude) vs latitude for 46-point N-S 
grid, wave number 8, phase speed 10°/day and A = 1.6 x 10 3 . 
(Latitudes with zero amplitude are not included and time is 
given in hours.) 
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Figure 22. Terrain pressure amplitude vs latitude for initial field 
and 48-hour forecast, wave number 8, phase speed 10° /day 
and A = 1.6 x 10 8 . 
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Chart I. 



Initial surface pressure analysis for 46-point N-S grid, 
wave number 8, phase speed 10°/day and A = 1.6 x 10 e . 
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Chart J. 48-hour surface pressure forecast for 46-point N-S grid, 
wave number 8, phase speed 10° /day and A = 1.6 * 10 8 . 
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PHASE ANG IE (DEG RE ES LONGITUDE) 

Figure 23. Phase angle (degrees longitude) vs latitude for 46-point N-S 
grid, wave number 12, phase speed 10° /day and A = 7.0 * 10 7 . 
(Latitudes with zero amplitude are not included and time is 
given in hours.) 
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Figure 24. Terrain pressure amplitude vs latitude for initial field 
and 48- hour forecast, wave number 12, phase speed 10° /day 
and A = 7.0 x 10 7 . 
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Chart K. Initial surface pressure analysis for 46-point N-S grid, 
wave number 12, phase speed 10° /day and A = 7.0 x 10 7 . 
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Chart L. 48-hour surface pressure forecast for 46-point N-S grid, 
wave number 12, phase speed 10° /day and A = 7.0 x 10 7 . 
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VI. CONCLUSION 



The model remained well-behaved throughout all experi- 
ments and the phase speeds were less than for the non-diver- 
gent model, however, they equaled or exceeded theoretical 
expectations as described by Williams (1972). The model, as 
it is presently constructed, appears to have the desired 
horizontal flexibility as evidenced by the fact that the ex- 
periments were carried out with different grid spacings. 

The vertical flexibility, although built in, was not tested 
in this report. The simplified nature of the model pre- 
cluded direct comparison with existing global models such as 
FNWC's 5-level model. 

Future research efforts with this model should most 
likely include expansion to 5-levels, treatment of cross 
polar flow and initialization with real data. 
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